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Abstract. We discuss the efficiency of Monte Carlo methods in solving continuum radiative transfer problems. The sampling 
of the radiation field and convergence of dust temperature calculations in the case of optically thick clouds are both studied. 
For spherically symmetric clouds we find that the computational cost of Monte Carlo simulations can be reduced, in some 
cases by orders of magnitude, with simple importance weighting schemes. This is particularly true for models consisting of 
cells of different sizes for which the run times would otherwise be determined by the size of the smallest cell. We present 
a new idea of extending importance weighting to scattered photons. This is found to be useful in calculations of scattered 
flux and could be important for three-dimensional models when observed intensity is needed only for one general direction 
of observations. Convergence of dust temperature calculations is studied for models with optical depths tv = 10 — 10 4 . We 
examine acceleration methods where radiative interactions inside a cell or between neighbouring cells are treated explicitly. In 
optically thick clouds with strong self-coupling between dust temperatures the run times can be reduced by more than one order 
of magnitude. The use of a reference field was also examined. This eliminates the need for repeating simulation of constant 
sources (e.g., background radiation) after the first iteration and significantly reduces sampling errors. The applicability of the 
methods for three-dimensional models is discussed. 
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1. Introduction 

Continuum radiative transfer problems are commonly solved 
with Monte Carlo calculations. The radiation field is simulated 
according to the actual processes that are assumed to be tak- 
ing place in an interstellar cloud. Emission from background, 
internal sources, and the dust itself is represented by a num- 
ber of photons packages, each corresponding to a large number 
of actual photons at one wavelength. Random numbers deter- 
mine the initial positions and directions of the photon pack- 
ages. The distance to a point where scattering takes place and 
the direction after scattering are also calculated using random 
numbers. This is the basic scheme in studies of light scattering 
(e.g., Mattila [1970; Witt 119771 . In scattering, the probability 
distribution for the change of photon direction is given by the 
scattering function that depends on dust properties. It is this 
random change of direction that precludes (with the exception 
of isotropic scattering or pure forward scattering) the use of 
direct ray tracing methods which are often used in line trans- 
fer. If, in addition to scattered light, also dust emission is to be 
solved then during the simulation the number of absorbed pho- 
tons must be registered in each cell of the model cloud. This 
could be done at each position where scattering takes place but 
unless optical depth is very high (or one is calculating only the 
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scattered flux) it is better to explicitly calculate absorptions in 
all cells that the photon package passes through. Especially in 
optically thin clouds the statistics of absorbed energy would 
otherwise remain very poor. 

Monte Carlo continues to be the standard method for the 
modelling of dust emission from interstellar clouds (Bernard et 
al. 1 19921 for recent papers see, e.g., Stamatellos & Whitworth 
120031 Concalv es et al. 120041 Niccolini et al l2003l Juvela & 
Padoan 2003 Pascucci et al. 120041 and references therein; 
Kurosawa et al. l2004l Whitney et al. l2003l Wolf et al. l2003l . 
Several methods can be employed to improve the sampling. 
One of the most common is the method of forced first scatter- 
ing (e.g., Mattila 1970 1 which improves sampling of scattered 
flux in clouds of low optical depth. In an optically thin cloud 
most photons pass through the cloud unscattered and give no 
information on scattered flux. Forced scattering means that one 
calculates for each package the fraction of photons that do scat- 
ter in the cloud. Unscattered photons are followed through the 
cloud along the original direction while the calculated fraction 
is always scattered somewhere along the way. The distance to 
the point of scattering is calculated from a conditional proba- 
bility distribution 'where will the photon package scatter if it 
does scatter before exiting the cloud' . 

In optically thick clouds most of the incoming flux is scat- 
tered many times and a photon package does not usually propa- 
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gate very far. In a spherically symmetric model this means that 
very few background photon packages reach the cloud centre. 
In order to get proper sampling in optically thick regions the 
number of simulated photon packages must be extremely large, 
and the run times become correspondingly very long. Niccolini 
et al. (120031 method 2) presented a partial solution which is 
analogous to the method of forced scattering but where the 
roles of the scattered and unscattered flux are reversed. The 
fraction of photons that does not scatter is first calculated and 
that part of the photon package is followed through the cloud 
toward the original direction. Rest of the photons is scattered 
somewhere along the way. The procedure is repeated after each 
scattering. While the method of forced scattering gives (in sta- 
tistical sense) the right amount of scattered flux it is not equally 
clear that the method of Niccolini et al. ( 2003 1 always provides 
correct estimates of absorbed photons. If the number of simu- 
lated photon packages is small the flux observed in cloud centre 
consist only of photons that are first scattered a few times near 
the surface and after this propagate unscattered to the centre. It 
remains open what is the ratio of these photons and those in the 
more rare photon packages (not necessarily more rare photons) 
that reach the centre only after many scatterings. At the limit of 
high optical depth for scattering and low optical depth for ab- 
sorption the method would appear to improve the sampling but 
actually all flux in the centre should still result from packages 
that have scattered numerous times, i.e. the method could not 
decrease the number of photon packages needed in the simu- 
lation. These are, of course, not faults of the method itself but 
are simply consequences of the sampling problem. The method 
does not address problems arising from different cell sizes and, 
therefore, while this method is probably very useful in many 
cases it is not yet a complete solution to the problem. 

One of the main advantages of Monte Carlo methods is 
their great flexibility in the way radiation field is sampled. All 
probability distributions employed in the simulation of photon 
packages can be modified as long as these are compensated 
by corresponding changes in the weighting, i.e. the number of 
actual photons in a package. For some reason this advantage 
is often not used and simulations are unnecessarily inefficient. 
High optical depths are perceived to pose another problem and 
not only because of difficulties in the sampling of the radia- 
tion field. In an optically thick cloud there may be significant 
self-coupling between dust temperatures. In the usual scheme 
this means that dust temperatures can not be solved directly 
and one is forced to iterate between simulation of the radiation 
field and updating of dust temperature estimates. The number 
of iterations depends on the optical depth, r. For a cloud where 
t reaches at optical wavelengths several hundred run times can 
become orders of magnitude longer than in the case of an op- 
tically thin cloud. This depends, however, very much on the 
external radiation field and the dust temperatures. According 
to Bernard et al. |1992J), for cold clouds embedded in normal 
interstellar radiation field the dust coupling is unimportant un- 
less optical depth is several hundred. This was confirmed by 
Concalves et al. (2004; see their Fig. 1) but the situation can be 
very different in the presence of hot dust. 

The slow convergence at high optical depths is, of course, 
not related to the Monte Carlo method per ce. Even if radiation 



field is sampled with Monte Carlo method the dust tempera- 
tures can be solved without iterations and even if iterations are 
made, simulation of the radiation field need not be repeated. 
Because of the larger memory requirements such schemes are, 
however, practical only for one-dimensional models. Bjorkman 
& Wood ( 2001 ) presented a modification where local dust tem- 
perature is updated after each absorption and a new photon is 
immediately re-emitted from the same position. The frequency 
of the emitted photon is obtained from a probability distribu- 
tion that takes into account both the previously emitted photons 
and the correction due to the updated temperature. The desired 
noise level determines the number of required photon pack- 
ages and once these and the induced re-emitted photons have 
been simulated the calculations provide a converged solution. 
The authors state that the method does not suffer from conver- 
gence problems associated with A-iteration methods. We return 
to this question later in the paper. 

In this paper we will study the efficiency of Monte Carlo 
radiative transfer calculation. We will discuss separately ef- 
ficient sampling of the radiation field and convergence of it- 
erations. One-dimensional, i.e. spherically symmetric model 
clouds are used as examples. We will, however, emphasise that 
most results can be directly transferred to three-dimensional 
cloud modelling and we will briefly discuss the applicability of 
the methods in the case of 3D models containing millions of 
cells. In Sect. |2] we describe the program and the implementa- 
tion of weighting, reference field and convergence acceleration 
schemes. Results from tests with one-dimensional models are 
presented in Sect.|3j and the results are discussed in Sect.|4] 

2. The program 

In this section we present some details of our implementation 
of the radiative transfer program. More specifically, we discuss 
the possibilities of improving the sampling of the radiation field 
so that a given accuracy of the results can be obtained with a 
smaller number of photons packages. On the other hand, we 
will discuss the convergence of dust temperature calculations 
in the case of optically thick clouds. This part of the discussion 
does apply equally to programs where some method other than 
Monte Carlo is used to estimate the radiation field. 

2. 1. Sampling of the radiation field 

One of the main advantages of Monte Carlo methods is that 
the sampling can be easily adapted according to the problem at 
hand. This possibility is, however, seldom used and, in order to 
ensure proper sampling in all parts of the model, simulations 
need an exceedingly large number of photon packages. This 
adaptation or 'weighting' is what in Monte Carlo integration 
would be called importance sampling. While it may be possible 
to create adaptive methods (analogous to stratified sampling; 
see Press et al. U992> we restrict our discussion to more simple 
schemes. 

In normal Monte Carlo simulation photon packages are cre- 
ated at random locations (within emitting medium or on the 
surfaces of separate radiation sources) and sent uniformly to- 
ward random directions. Each package has equal weight i.e. 
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Fig. 1. The four weighting schemes (A, B, C, and D) used in 
the Monte Carlo simulations in this paper. The figures illustrate 
a case where the sampling is improved in the centre of a spher- 
ically symmetric model cloud: background packages are sent 
preferentially towards the cloud centre (scheme A), more pho- 
ton packages are created in the inner parts (scheme B), within 
the cloud photon packages are sent preferentially towards the 
cloud centre (scheme C), and scattered photon packages are 
directed preferentially towards the cloud centre (scheme D). 

the number of true photons is divided evenly between created 
packages. One is not restricted to the use of uniform spatial and 
angular distributions and one could, e.g., send twice as many 
photon packages from certain part of the cloud provided that 
each of these packages contains only half of the original num- 
ber of photons. More generally, if one assumes a new probabil- 
ity distribution p(r, ft) for positions r and directions of the 
created photon packages the weight of the photon package is 

' p(ri,Sli)' 

The default number of photons included in a package is mul- 
tiplied with this number where and fij are the actual po- 
sition and direction of the photon package. Weighting can be 
applied to background photons, photons emitted by dust in- 
side the cloud and for each discrete source included in the 
model. In each case the spatial and angular distributions can 
be weighted separately. The four weighting schemes used in 
this paper (A, B, C, and D) are illustrated in Fig.^ These are 
are discussed below and a more detailed presentation can be 
found in Appendix lA.21 

2.1 .1 . External radiation field: weighting scheme A 

In a spherically symmetric cloud all background photon pack- 
ages enter the outermost shell but, because of their smaller di- 



ameter, the inner shells are hit by relatively few packages. At 
least in optically thin clouds the problem can be alleviated by 
sending photon packages preferentially toward cloud centre. 
The original probability distribution of the impact parameter, 
p{d) oc d, can be replaced one with higher probability at small 
values of d (e.g., p(d) oc d a with a < 1). We use a scheme 
where the same number of photon packages is sent towards 
each annulus as defined by the radial discretization (weighting 
scheme A). Within each annulus we have the usual distribu- 
tion of p(d) oc d. The number of photon packages is fixed for 
each annulus and the random noise is reduced. This is, how- 
ever, not yet optimal. In an optically thin model with N shells 
each photon package enters the outermost shell but only one in 
N hits the innermost cell. A more strongly peaked probability 
distribution could result in more uniform errors. 

2.1.2. Internal radiation field: weighting scheme B 

Different cell sizes also affect the sampling of dust emission 
within the cloud. This is particularly problematic in optically 
thick clouds where self-coupling between dust temperatures 
is strong. In a spherical model the the innermost cells can 
be orders of magnitude smaller than the the outermost shells. 
Sampling errors increase toward centre and the number of re- 
quired photon packages depends directly on the size of the 
smallest cells. The solution is similar as in the case of back- 
ground photons. Emissions take place at random positions with 
radial probability distribution, p(r) ~ r 2 . For photon packages 
this can be replaced with a steeper distribution (e.g., p(r) ~ r@ 
with below 2). We use a procedure where the same number 
of photon packages is generated in each cell but the locations 
within the cells follow the p(r) ~ r 2 distribution (weighting 
scheme B). In extremely optically thick clouds one might also 
consider whether photon packages should be created preferen- 
tially close to cell boundaries so that these would better sample 
the energy transfer between cells. The importance of different 
cells as sources of radiation could also be taken into account 
and more packages could be sent from denser and warmer cells. 

2.1 .3. Angular distribution of photons created in the 
cloud: weighting scheme C 

In order to improve the sampling of the smallest cells we con- 
sider weighting the angular distribution of photon packages 
created within the cloud. Uniform angular distribution is re- 
placed with distribution that peaks in the direction of the cloud 
centre or, more generally, toward any region where sampling is 
to be improved. We use an exponential function p{9) oc e ,e 
where 8 is angular offset from the direction of the cloud cen- 
tre and 7 is a positive constant (weighting scheme C). The 
exponential function is convenient since the inverse cumula- 
tive probability density function, P^ 1 , is easily calculated and 
can be used to generate angles from the selected distribution. 
However, the calculation of the P^ 1 can be replaced with a 
lookup to a table of P-values and, in practise, there are no re- 
strictions for using any function as the probability distribution 
P- 



4 



M. Juvela: Efficient Monte Carlo methods 



2.1 .4. Place and angular distribution of scattering: 
weighting scheme D 

In weighting scheme D we apply weighting to scattered pho- 
tons. The method of forced first scattering would be one exam- 
ple of this (see Sect. IA.U . Here we are, however, more inter- 
ested in the angular distribution of the scattered photons. The 
goal can be to increase the flux of packages (not photons!) to- 
ward optically thick parts of the cloud or toward one partic- 
ular direction from which the emerging scattered intensity is 
observed. Scattered photons already have a non-uniform prob- 
ability distribution, p(SY), which is determined by scattering 
function of the dust model. This can be given by the Henyey- 
Greenstein dl941> formula, or some other function specified 
in analytical or tabular form. The distribution p is defined in 
a coordinate system that is fixed by the initial direction of the 
photon package. Let us assume that we want to have locally 
a probability distribution q(fl) for the directions of scattered 
photon packages. Random directions are generated according 
to this distribution and each package is weighted with the ratio 

W, - *§. (2) 

Functions p and q are defined in different coordinate systems 
and one must apply a rotation to find out what direction fl' the 
selected direction £1 corresponds to. 

The presented scheme clearly works if one divides a pho- 
ton package into n small packages, each containing a fraction 
oc p(O')/ q(fi) of the original number of photons and having di- 
rections following the distribution q. For practical reasons it is 
not possible to split a package this way. One can, however, use 
the scheme without creating any new photon packages. When 
a package is scattered the new direction is obtained from the 
probability distribution q and the number of photons contained 
in the package is multiplied by the factor given by Eq. [2] The 
weight Wi can be smaller or larger than one and for some pho- 
ton packages the number of photons can actually increase. The 
result is, of course, correct only statistically and after the sim- 
ulation of many photon packages. The possibility of weighting 
leads to an interesting conclusion that all frequencies can be 
calculated simultaneously using a single photon package. Each 
time scattering occurs the number of photons at each frequency 
is multiplied by a factor depending on the realised scattering 
angle and the scattering function at that frequency. Weighting 
should also be applied according to the distance between scat- 
terings which depends on the frequency-dependent opacity of 
the medium. Since this is very wavelength dependent such a 
scheme would, however, most likely result in increased sam- 
pling errors. 

The method D is not the same as the 'peeling off method 
of Yusef-Zadeh et al. ( 1994 1 that is used to improve the quality 
of images of scattered light (Wood & Reynolds 2003 1 Whitney 
et al. 2003 1. In the 'peeling off method one calculates and reg- 
isters after each scattering the fraction of photons that are scat- 
tered towards the observed and escape the cloud. The angular 
distribution (or weighting) of the rest of the photons still essen- 
tially follows the original scattering function. Therefore, the 
sampling of the scattered flux remains unchanged within the 



cloud while in our method D it is changed. Method D improves 
the sampling of the observed scattered flux by increasing the 
number of regular photon packages that exit the cloud towards 
the observer. The 'peeling off method is more efficient if only 
emerging flux is considered and only for one single direction. 
Method D might be competitive in cases where observed flux 
results from multiple scatterings in an optically thick medium. 
Because method D improves sampling in the selected general 
direction it can also be used when images of the scattered flux 
are needed for several directions close to each other. Of course, 
it would be possible to combine the two methods. Method D 
would be used to modify sampling within the cloud and the 
'peeling off could still be performed after each scattering. This 
possibility is not studied any further in this paper. 

2.2. Reference field 

In Monte Carlo integration errors are proportional to function 
values. If solution is known for a reference function that is sim- 
ilar to the integrated function only the difference needs to be 
solved with Monte Carlo methods and the errors are propor- 
tional to the difference between the actual function and the ref- 
erence. This idea was first applied to line transfer calculations 
by Bernes (1979). He assumed a reference field corresponding 
to a fixed excitation temperature and Monte Carlo simulation 
was used to determine the differences between the true field 
and the constant reference field. Choi et al. (1994) improved 
the method by using the solution from the previous iteration as 
the reference. As the reference field approaches the true field 
the sampling errors decrease on each iteration. 

So far a reference field has not been used in connection 
with dust emission calculations. The implementation is, how- 
ever, straightforward. The reference field is taken to correspond 
to the situation of the previous iteration. On the first iteration 
reference field is zero and calculations proceed in the normal 
fashion. On the following iterations each created photon pack- 
age contains the difference between the true number of pho- 
tons (corresponding to the latest temperature estimates) and the 
photons from the reference field (corresponding to the previous 
temperature estimates). For constant radiation sources and the 
background radiation this difference is zero and simulation is 
needed only to find out the effect of the latest dust temperature 
updates. 

The concept of a reference field is useful in calculations 
involving many iterations. First of all, emission from constant 
sources (e.g., background) can be simulated on a single itera- 
tion. On the following iterations their effect would already be 
included in the reference field and no further photon packages 
need to be simulated from them. Secondly, sampling errors are 
decreased (or a given accuracy is reached with fewer photon 
packages) since the final noise level depends on the total num- 
ber of simulated photon packages rather than the number of 
packages per one iteration. Details of the implementation are 
given in Appendix lA.41 In principle, the use of a reference field 
does not affect the convergence of dust temperatures that is dis- 
cussed below. 
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2.3. Accelerated iterations 

Only in the optically thin case (t << 1) can dust tempera- 
tures be solved directly. In optically thick clouds the dust emis- 
sion can contribute significantly to dust heating and this leads 
to iterations where radiation field and dust temperatures are 
solved alternatingly. As r increases, more and more iterations 
are needed and calculations can potentially become very time- 
consuming. 

Dust temperatures can, in principle, be computed without 
repeating the simulation of the radiation field. Once the ra- 
diative coupling between all cells has been calculated one can 
write a set of equations from which dust temperature in all cells 
can be solved. The same applies to, e.g., molecular line calcu- 
lations where the unknowns are level populations in different 
cells. As noted by Stamatellos et al. (2003 1, the case is much 
simpler for the continuum radiative transfer in dust clouds. 
Absorption and scattering cross sections can be assumed to be 
independent of temperature and radiative coupling (e.g., 'the 
fraction of photons emitted from cell i is absorbed in cell f) 
remains constant. If information about the coupling is saved 
a simulation providing this information needs to be done only 
once. If one considers dust particles at equilibrium tempera- 
ture the whole problem reduces to a set of non-linear equations. 
The problem becomes more complicated if one must consider 
in each cell a distribution of dust temperatures. Temperature 
distribution is usually solved separately from a linear set of 
equations where the unknowns (i.e. number of dust particles 
in a given temperature interval) are multiplied with factors that 
depend on temperature distributions in other cells. The whole 
problem (including radiative couplings) can still be formulated 
as a single set of non-linear equations. 

For one-dimensional problems where only dust particles at 
an equilibrium temperature are considered the solution of such 
equations is still feasible. For transiently heated particles or for 
any models in several dimensions this is no longer practical. 
For a model of N cells one would have to store for each sim- 
ulated frequency (N 2 — N)/2 terms describing the coupling 
between any two cells. If these can not be stored then also the 
simulation of the radiation field needs to be repeated each time 
before temperatures can be updated. For optically thin clouds 
iterations are never needed. Once the optical depth becomes 
high hundreds of iterations may be required and it becomes 
difficult to follow the convergence if results include random 
sampling errors. The slowness of convergence is, of course, not 
caused by the Monte Carlo simulation which is only a method 
for estimating the radiation field. 

As r increases most emitted photons are absorbed locally 
and a larger fraction of energy flow takes place within a cell 
with relatively small amount of interaction between cells. In 
calculations this translates to small temperature corrections and 
slow convergence of dust temperatures. The number of itera- 
tions needed to reach correct dust temperatures depends not 
only on optical depth but also on dust temperatures and the 
spectrum of the external radiation field. If dust is cold and, in 
spite of high visual optical depth, cloud is penetrated with suf- 
ficient amount of longer wavelength external radiation the self 
coupling of dust temperatures may remain unimportant. 



There are several ways to improve the convergence rate. We 
will consider two alternatives, a heuristic extrapolation based 
on dust temperatures adopted from previous iterations, and the 
accelerated Monte Carlo methods. Accelerated Monte Carlo 
schemes (AMC) are analogous to ALI (Accelerated Lambda 
Iteration) methods and in these part of radiative interactions 
are treated explicitly when dust temperatures are updated. Such 
methods have already been used in line transfer (Juvela & 
Padoan l200Tl Hogerheijde & van der Tak:2000). We will con- 
centrate on the case where we can assume that dust is in each 
cell at one equilibrium temperature. The generalisation to the 
case of a dust temperature distribution is rather straightforward. 

2.3.1 . Heuristic extrapolation 

Extrapolation is based on dust temperatures at three con- 
secutive iterations. We denote the latest change per iteration 
(AT/iter)i, and the factor k by which this derivative was seen 
to change, i.e. (AT /iter), = k(AT/iter) l _ 1 . If < |fc| < 1, 
i.e. the temperatures show signs of convergence, an extrapo- 
lation is done. This assumes that each successive temperature 
correction will continue to be k times the previous one. The 
extrapolation approaches a constant value but in our imple- 
mentation we make extrapolation only 20 iterations forward. 
The temperature correction obtained from the extrapolation is 
further limited to maximum of 20% of the original value. The 
method requires some extra memory (two values per cell) but 
very little computations. In the simulation the same set of ran- 
dom numbers should be used on each iteration so that random 
variations do not interfere with the extrapolation. 

2.3.2. AMC methods 

In normal calculations one must compute in each cell the 
strength of the radiation field which consists of the intensity 
produced by the cell itself and the intensity caused by every 
other cell. Later this information is used to compute new esti- 
mates of dust temperatures. The other extreme case was men- 
tioned above: the coupling between cells is determined and 
temperatures are computed simultaneously from a large set of 
equations. In between there are many alternatives where part of 
the interactions are treated explicitly in dust temperature cal- 
culations. This is analogous to Accelerated Lambda Methods 
(ALI) used in line transfer (Rybicki & Hummer [T99T1IT9921 . 
There, the local intensity J is computed from source functions 
S through lambda operator, J = AS, and the operator is split 
into two parts. In the most common case the diagonal part of 
the operator is separated. This represents the intensity produced 
in a cell by the cell itself. In an optically thick case most of the 
created photons are absorbed locally. This loop does not con- 
tribute to a change in level populations but its elimination does 
significantly improve the rate of convergence. 

In Monte Carlo we calculate photon absorptions rather than 
the intensity. The use of diagonal A operator is replicated by 
calculating in each cell only those absorbed photons that were 
emitted outside that cell. At the same time we note for each 
cell the photon escape probability i.e. the fraction of photons 
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that escapes the emitting cell. One must remember to take into 
account those photons that first escape the cell but are later 
scattered back and absorbed. In dust temperature calculations 
absorbed photons originating from other cells, discrete radi- 
ation sources and from the background are balanced against 
that fraction of emitted photons that leave the emitting cell. 
Temperature calculations do not become more complicated but 
some additional memory is required as the photon escape prob- 
abilities must be saved for each cell and each simulated fre- 
quency. In spite of this the method is feasible even for three- 
dimensional models. We will call this the AMC-method with 
diagonal operator. 

We have also implemented a 'tridiagonal operator'. For 
one-dimensional models this means that when dust tempera- 
tures are updated the radiative coupling between closest neigh- 
bours is treated explicitly. In simulations we ignore those ab- 
sorbed photons that were emitted by the cell itself or one of 
its immediate neighbours. Additionally, we note again the ra- 
diative coupling between neighbouring cells i.e. the fraction 
of emitted photons that is absorbed in each of the neighbour- 
ing cells. Dust temperatures are solved from a non-linear set 
of equations where the i:th equation contains three unknown 
temperatures, one for the i:th cell and one for each of its neigh- 
bours. We have used a simple iteration where one temperature 
at a time is updated until all temperatures have converged. The 
calculation converges with very few iterations, and time spent 
in solving these equations is insignificant compared with the 
overall run-times. The use of tridiagonal operators is not nec- 
essarily restricted to one-dimensional models. It could be used 
in three-dimensional models where, due to the geometry, each 
cell interacts mainly with two neighbours (e.g., pieces of spher- 
ical shells). 

The implementation of the AMC methods is described in 
more detail in Appendix lA.3l 

3. Tests with one-dimensional models 

In this section we test in practise the methods listed in the pre- 
vious section. Improvements in the sampling of the radiation 
field and in the convergence of dust temperature calculations 
are studied separately. Tests are made using spherically sym- 
metric model clouds and assuming dust at an equilibrium tem- 
perature. 

The first model is taken from Stamatellos et al. ( 2003 their 
model BE2). The cloud represents a sub-critical, externally 
heated spherical cloud. The radial density distribution follows 
the Bonnor-Ebert solution for an isothermal spherical cloud 
bounded by external pressure. In connection with this cloud we 
use the dust model of Ossenkopf & Henning ( 1994; coagulated 
grains with thin ice mantles accreted in 10 5 years at a density 
of 10 6 cm~ 3 ) and the external radiation field as given by Black 
( 1994 1. Therefore, our calculations correspond exactly to those 
presented by Stamatellos et al. The visual optical depth to the 
centre of the cloud is 16.6. In the following this is called the 
model Si. We study also an optically thicker cloud, model S*2, 
where the optical depth has simply been scaled up by a factor 
of ten. Model clouds are divided into 50 shells of equal geo- 



metric thickness and the number of simulated frequencies was 
44. 

From Niccolini et al. (2003 1 we take a model where the 
optical depth to cloud centre is 10 at 1 /im. The cloud has a 
central radiation source with a spectrum corresponding to a 
black body with T=2500K. We study also a second model, N2, 
which is identical to Ni but has an optical depth r(l ^m)=100. 
The models consist of 51 shells. The radiae are equidistant on 
logarithmic scale except for some very narrow shells close to 
the cavity surrounding the central source. In accordance with 
Niccolini et al. (2003 1 we use a simple dust model where ab- 
sorption and scattering cross sections are proportional to the 
frequency and scattering is isotropic. The number of simulated 
frequencies was 40. 

3. 1. Sampling in temperature calculations 

Fig-EJshows convergence as the function of the number of sim- 
ulated photon packages per iteration (and frequency), n p , for 
model Si . Maximum error and overall rms error are shown for 
derived dust temperatures in different shells. The errors were 
determined by comparison with a calculation with much higher 
rip. Results are plotted for normal Monte Carlo sampling (no 
weighting) and for calculations where equal numbers of exter- 
nal photon packages were sent towards each annulus described 
by the radial discretization (weighting scheme A) and equal 
number of photons packages were sent from within each cell 
(scheme B). Equal number of photon packages were used for 
describing the external field and emission within the cloud. The 
convergence is for both methods ~ 1/ ^/n^ and differences in 
accuracy remain relatively small. However, for weighted sam- 
pling the errors are always equal or smaller than in normal runs. 
In particular, fluctuations of the results are smaller and for some 
n p the rms error is almost one order of magnitude smaller than 
in runs where weighting was not applied. 

Fig. [3] shows same relations for model S2 that has a fac- 
tor of ten higher optical depth. In this case the advantage of 
weighted sampling is more noticeable. On the average, errors 
have decreased 'only' by a factor of ^8 but this translates to an 
almost two orders of magnitude difference in run times. 

Models Ni and N2 are fundamentally different because of 
the internal source. Both models are optically thick and dust 
temperatures are mostly determined by dust emission itself. 
For optically thin models the default scheme could be improved 
only by making sure that photon packages are sent as uniformly 
as possible towards different directions. This could be accom- 
plished with the use of quasi random numbers or equally by 
using pre-selected, angularly equidistant directions. 

The optical depth towards the centre of cloud Ni is 10 at 
1/U and the visual extinction is roughly twice this value. The 
angular distribution of photons emitted from the source is not 
weighted, background photons are not included at all and we 
test only the effect of weighting the distribution of locations 
at which dust emitted photon packages are created within the 
cloud (method E). Results are shown in Fig.|4]where relative 
rms-errors from un-weighted and weighted runs are compared. 
If photon packages are created uniformly over the cloud vol- 
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Fig. 2. Convergence of calculated dust temperatures in model 
Si as the function of the number of photon packages per 
iteration and frequency. Results are shown for both normal 
Monte Carlo calculations (filled squares) and calculations 
where weighting was applied to both dust emission and back- 
ground emission (open circles; see text for details). The upper 
curves show the maximum relative error in any of the shells 
and the lower curves show the rms-value of the relative errors 
summed over all shells. The dotted line shows the expected 
convergence rate. 
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Fig. 3. Convergence of calculated dust temperatures in model 
S*2 vs. number of photon packages per iteration and frequency. 
The symbols are as in Fig.|2] 



ume the sampling of the innermost shells remains very poor. 
The results are essentially incorrect unless the number of pho- 
ton packages is larger than the ratio between the volume of the 
cloud and the volume of the smallest cells or rather the vol- 
ume of smallest optically thick region that could separate the 
outer cloud from the source. The observed convergence of un- 
weighted results indicates that the thickness of this inner re- 
gion is in this case a few per cent so that only about one photon 
package out of 10 5 samples emission from there. Results show 
clearly the need of some kind of weighted sampling. 

In model iV2 optical depths are higher by a factor of ten 
and normal Monte Carlo sampling becomes even more inef- 
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Fig. 4. Rms-value of relative dust temperature errors as func- 
tion of the number of simulated photon packages per iteration 
and frequency. Results are shown for model TYi for the nor- 
mal Monte Carlo simulations (filled squares) and for weighting 
scheme B where relatively more photon packages are created 
close to the centre of the model (open circles; see text for de- 
tails). 

ficient. This is clearly seen in Fig. |5} where no convergence 
is seen before the number of photon packages is ~10 6 and 
adequate sampling of the inner region becomes possible. For 
the weighted sampling the situation is much better and actu- 
ally not worse than in the case of the previous model N\. As 
a result, a given accuracy is achieved with a number of photon 
packages that is a factor of 10 4 lower than in the case of un- 
weighted sampling. Another remarkable fact is that the conver- 
gence seems to be much faster than the usual 1 / v^V behaviour. 
This is not altogether surprising since for regular sampling (as 
accomplished in Monte Carlo integration by the use of quasi- 
random numbers; see eg. Press & Teukolsky 19921 Il989b the 
convergence is expected to be ~ 1/N rather than ~ 1/^/N. 
In our case the regularity is restricted to systematic selection 
of shells and while positions and directions within a shell are 
completely random this is enough to ensure the faster conver- 
gence rate. 

Previously, we used two weighting schemes ( A and B) in 
comparison with the default Monte Carlo sampling. As dis- 
cussed in Sect. I2.1l there are at least two further ways to influ- 
ence the sampling, especially in the cloud centre where there 
would otherwise be very few photon packages. One can send 
photon packages from cells preferentially towards the cloud 
centre (scheme C), and scattered photon packages can be di- 
rected preferentially in that direction (scheme D). We tested 
the effect of these methods for models Si, S2, N\, and N2, 
with the number of photons packages (per iteration and fre- 
quency) equal to n p = 2000. The rms values of the relative 
temperature errors over all cells are shown in Table[0for some 
combinations of the four weighting schemes. 

For the externally heated models Si and S2 the scheme A 
(weighted generation of background photon packages) was the 
most useful one and it reduced rms errors by a factor of ^3 in 
model Si and by a factor of ~5 in model S2. Methods B and 
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Fig. 5. Convergence of dust temperature calculations in model 
N2 as the function of the number of photon packages per itera- 
tion and frequency (symbols as in Fig. 

Table 1. Errors of computed temperature for models S\, S2, 
N\, and N2 when default Monte Carlo sampling or different 
weighting schemes (A, B, C, and D) are used. The numbers 
are rms values of relative errors given as percentages. Number 
of photon packages per iteration was 2000. 



weighting scheme 


Si 
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Ni 


N 2 


none 


3.2 


1.7 


34 


81 
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1.1 


0.34 
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3.2 


1.7 


0.49 


4.7 
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3.4 


1.4 


36 


81 
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1.7 


1.7 


37 


82 


A,B 


1.1 


0.46 






A,B,C 


1.1 


0.31 






A,B,D 


0.46 


0.48 






A,B,C,D 


0.46 


0.31 







C (distribution of positions and directions of photon packages 
created inside the cloud) affected the results only little. This is 
not surprising since in these models the dust emission has but 
a small contribution to the radiation field. Method C produced 
small improvements only in cloud S2 while method D was ef- 
fective only in model Si. The best combinations were A, B, 
and D in model Si and A, B, and C in model S2, and the rms 
errors were correspondingly reduced by factors ~ 7 and ~ 6. 

In internally heated clouds Ni and N2 the emission from 
innermost small and hot cells must be properly sampled - oth- 
erwise the flow of the re-emitted energy through the optically 
thick cloud is cut and results will be incorrect. Table ^ shows 
that errors are very large if scheme B is not used. Scheme A 
was not applied since background photons were not simulated. 
The other weighting methods have no real effect on the errors. 
Methods C and D aim at improving in cloud centre the sam- 
pling of radiation that originates further out. In clouds Ni and 
N2 the flux of energy is mostly in the reversed direction and 
the use of these methods is not helpful. 



3.2. Sampling of scattered photons 

The weighting of the angular distribution of scattered photon 
packages (scheme D) can be used to improve sampling in se- 
lected regions inside the cloud. It can also be used to improve 
the statistics of scattered flux observed outside the cloud. For 
three-dimensional models the out-coming flux might be needed 
only for one direction. In that case, scheme D can be used to 
drive photon packages towards the selected direction. The sta- 
tistical noise will be decreased for the scattered flux and the 
same simulation can still be used for solving the temperature 
structure of the cloud. 

We calculated for model S2 400 pixel maps of scattered 
flux at 0.55/im. The method of forced first scattering was used 
in all calculations. In the weighting scheme D the probability 
function of the angular distribution of scattered photons was 
p ~ exp(—j6) where the angle 6 is measured from a vector 
pointing towards the observer. A reference solution was ob- 
tained with normal Monte Carlo using 4x 10 7 photons. In this 
map the expected error per pixel is ~2%. Runs with 2xl0 5 
photon packages were compared with this. The normal Monte 
Carlo method resulted in a standard deviation of 34%. The cor- 
responding error when using scheme D was ~ 20% when the 
parameter 7 was in the range 0.3-1.0. For a given noise require- 
ment and 1 / \/N noise dependence the difference corresponds 
to almost a factor of three in run times. One can expect that 
this kind of weighting could be more useful in the case of non- 
isotropic scattering or a non-isotropic radiation field. 

3.3. Convergence tests 

The methods of Sect 12.31 and their effect on the convergence 
of dust temperature calculations were tested first with model 
iVjj. In that model the optical depth to the cloud centre is 
t(1 /xm)=100. Calculations were started with a flat tempera- 
ture profile, T — 17K, while in the final solution the tem- 
peratures range from ~ 140 K to over 2000 K. Random number 
generators were reset after each iteration and the results were 
compared with similar calculations with a very large number 
of iterations. Therefore, the results do not show any sampling 
errors and reflect only the true convergence of dust temperature 
values. 

Fig-EDshows the rms-errors as function of the number of it- 
erations for two runs. With the normal method the convergence 
is quite slow and a relative rms error of 1 % is reached only after 
~600 iterations. When the extrapolation method of Sect. 12. 3.11 
is used the convergence improves by a factor of three and an 
accuracy of 1% is reached after some 200 iterations. The ex- 
trapolation step was done always after three normal iterations. 
The convergence is quite smooth and extrapolations did not at 
any point produce noticeable oscillations. This suggests that 
the convergence rate could be still improved by amplifying the 
computed corrections, at least during the first iterations. 

Results from corresponding tests with the AMC-methods 
are shown in Fig.0 With the diagonal operator (i.e. when in- 
ternal absorptions are treated explicitly) the convergence rate is 
about the same as with the extrapolation method and a relative 
accuracy of 1% is reached with ~200 iterations. When the two 
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Fig. 6. Convergence of dust temperatures in model iVjj in the 
case of normal iterations and with the use of the extrapolation 
method. The plotted values are rms-values of the relative tem- 
perature errors. 

methods are applied simultaneously this accuracy is reached 
with 60 iterations. The convergence rate is very good both ini- 
tially and again after some 45 iterations. In between the rate 
is slower, possibly because of some less successful extrapola- 
tion steps. In this test the tridiagonal operator is clearly better 
with a good and constant convergence rate. Finally, the circles 
in Fig. show the rms error after every fifth iteration when, 
in addition to the tridiagonal operator, an extrapolation step is 
done after every fifth iteration. This gives the fastest conver- 
gence and an rms error of 1% is reached with just 15 iterations, 
a factor of 40 less than for the basic method. 

3.4. Reference field 

The reference field concept was tested with model Ni which is 
more optically thin than N 2 but where the solution still requires 
many iterations. Figure [8] shows the results obtained with and 
without accelerated Monte Carlo and with and without a refer- 
ence field. In all calculations the number of photon packages 
per iteration and frequency was 4080, half of which described 
flux from the central source while the other half described the 
dust emission. The reference solution corresponded this time to 
calculations with both a larger number of iterations and photon 
packages per iteration. Because of the smaller optical depth the 
advantage of accelerated Monte Carlo is smaller than in model 
A<2 (see Fig. [7] above) but still clear. If the reference field is 
not used the final accuracy depends on sampling errors i.e. on 
the number of photon packages per iteration. When a reference 
field was not employed the same set of random numbers was 
used on each iteration. This explains the smoothness of the hor- 
izontal curve. If different random numbers were used the error 
level would vary from iteration to iteration but the general level 
would remain the same. 



Iterations 

Fig. 7. Convergence of dust temperatures in model N 2 for di- 
agonal AMC-method (solid line), diagonal AMC with extrap- 
olation (dotted line) and AMC with tridiagonal operator (dash- 
dotted line). 
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Fig. 8. Convergence of dust temperatures in model Ni with 
(filled symbols) and without (open symbols) the use of a refer- 
ence field when the number of photon packages per iteration is 
equal. The rms-values of relative temperature errors are plotted 
as function of the number of iterations. 

In principle, the use of a reference field does not affect the 
rate of convergence. However, after the first iteration no emis- 
sion from the central source is simulated since that information 
is already contained in the reference field. The following itera- 
tions should be about twice as fast if the overhead from the use 
of a reference field is ignored. In fact, measured after ten iter- 
ations the ratio of run times is 2.9:1 in favour of the reference 
field. Photon packages emitted from the central source are in- 
side the most optically thick region and therefore take a longer 
time to simulate. This explains why the ratio is in this case even 
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larger than 2:1. The second advantage of a reference field is that 
the sampling improves from iteration to iteration and the final 
errors depend on the total number photon packages simulated. 
In the case of Fig.|8]this means that convergence remains linear 
below the level that was reached without the use of a reference 
field. Emission from the central source was in this case simu- 
lated with 2040 photon packages per frequency while for dust 
emission the number of packages was 2040 times the number 
of iterations. It is clear that in calculations without a reference 
field the accuracy was limited by the sampling of the dust emis- 
sion. 

When the reference field is used the improvement as a func- 
tion of the number of iterations depends both on convergence of 
temperatures and improvement of sampling. In the test shown 
in Fig.|8]the number of packages per iteration was clearly suf- 
ficient so that the improvement of accuracy did not at any point 
decelerate because of sampling errors. Similarly, 2040 photon 
packages from the central source at each frequency were suffi- 
cient so that those sampling errors are not yet visible in Fig. [8] 
One could easily use more photon packages to simulate emis- 
sion from sources and background since that simulation is done 
only once. If the number of photon packages describing dust 
emission is too low it slows down the convergence but does not 
necessarily prevent it. In order to minimise the run times one 
should use as small a number of photon packages as possible. 
In the case of model Ni the number could be brought down to 
^100 without significantly affecting the convergence rate seen 
in Fig. El 

4. Discussion 

4. 1. Monte Carlo sampling 

Several weighting schemes were used to improve the sampling 
so that a given accuracy could be reached with fewer photon 
packages and shorter run times. All four schemes A, B, C, 
and D were useful in some conditions. For optically thin or 
moderately thick and externally heated clouds the scheme A 
(weighting of spatial distribution of background photons) was 
the most useful one while for optically thicker and internally 
heated clouds the scheme B was by far the most important one. 
The overall reduction in the number of photon packages ranged 
from a factor of a few to over ^10 4 . The results are extremely 
model dependent. However, a brief analysis of the model struc- 
ture and the expected flow of energy within the model will give 
a good idea on what kind of a weighting scheme would be most 
useful. 

Stamatellos & Whitworth ( 2003 ) sent external photons into 
the model cloud from one location at the cloud edge with prob- 
ability distribution p(0) ~ sinOcosO. This distribution does 
not give a particularly good sampling towards the cloud cen- 
tre (6 w 7r). The weighting scheme A ensured that the effect 
of the external field could be calculated accurately in the cloud 
centre no matter how the radial discretization was done. This 
would become more important if the number of shells were in- 
creased or the size of the innermost shells otherwise decreased. 
The systematic sampling used in the method (equal number of 
packages towards each annulus) also decreases Monte Carlo 



noise. As a consequence, the run times become rather short. For 
Si a relative 1% accuracy of temperatures was reached with 
3000 photons per frequency. Half of the photon packages was 
sent from the background and half from cells within the model 
cloud. We did not try to optimise this ratio. With an 600MHz 
PHI computer the run time was ^8 seconds per iteration. 

The method A should be most important for optically thin 
clouds but it was also found useful for the optically thick cloud 
52. One explanation is that while short wavelength radiation 
is absorbed in the outer layers the cloud centre is heated by 
longer wavelength radiation for which the optical depth may 
be close to unity. The spatial distribution created for external 
photons at the surface of the cloud is not completely destroyed 
by scatterings and will still result in more accurate estimates of 
the heating in inner regions. Details depend, of course, on the 
dust model and the spectrum of the external radiation field. 

For internally heated models Ni and N2 the method B 
was indispensable. In the case of model N\ an accuracy of 1 % 
was reached with just a few hundred photon packages per fre- 
quency. With 10 3 packages per iteration and frequency the run 
time with a 600MHz PHI computer was down to ~8 seconds 
per iteration. The required number of iterations can be read 
from Fig. [8] If the reference field is used the run time drops 
by more than half and can be shortened further by decreas- 
ing the number of simulated photon packages. We repeated the 
calculations with accelerated Monte Carlo (diagonal operator) 
and using a reference field, with 100 photon packages per fre- 
quency per iteration for the dust emission and 1000 photon 
packages for the central source in the first iteration. An rms- 
error of 1% was reached after 10 iterations and the whole cal- 
culation took less than 12 seconds. If these run times are com- 
pared with Niccolini et al. 2003 one should also note that while 
in our runs the cloud was divided into 5 1 shells they used 65 
shells (and 20 additional shells in the empty inner cavity). 

Method C (weighting of angular distribution of dust emis- 
sion) did produce small improvements only in model 52 but 
method D was not particularly useful in any of the four mod- 
els. This is partially due to the selection of the test cases. For 
example, in models Ni and N2 the source is situated in the 
cloud centre and the flow of energy is outwards while method 
D tried to improve the sampling of the scattered flux flowing 
in the opposite direction. Method D could still be useful, e.g., 
in cases where a cloud is heated from outside by a nearby star. 

A change in the sampling can affect the average time it 
takes to simulate a photon package. In an optically thin and 
spherically symmetric model those photons that are sent to- 
wards cloud centre pass through most cells and take more time 
to simulate. Frequent scatterings make simulations more time 
consuming in optically thick parts of the clouds. Therefore, 
weighting schemes of Sect. 12.11 should be compared not only 
against the number of photon packages required but also 
against the actual run times. Some tests were made with the op- 
tically thick model N\. With weighting schemes A and B the 
average run times per package were only ^3% longer than in 
the case of unweighted Monte Carlo. The increase is insignif- 
icant compared with the reduction in the required number of 
photon packages. For model 51 the difference in run times be- 
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Fig. 9. Convergence as function of the number of photon pack- 
ages for model Si with normal (squares) and weighted Monte 
Carlo (circles). Results obtained with pseudo random numbers 
(open symbols) and quasi random numbers (filled symbols) are 
shown 



tween normal Monte Carlo and weighting scheme A were sim- 
ilarly about 3%. 

The use of quasi random numbers was shortly tested but 
no significant improvement was observed when compared with 
results obtained with pseudo random numbers. For the un- 
weighted case the results from separate runs became more con- 
sistent (i.e. fluctuations decreased as seen in, e.g., Fig.|2}. The 
average error level was not improved and convergence rate (as 
a function of the number of photon packages) remained propor- 
tional to 1/ y 7 ?!^. The use of quasi random numbers is straight- 
forward for photon packages coming from the background and 
from internal sources. For emission within the cloud the posi- 
tion and direction of the created photon packages must be coor- 
dinated i.e. these should be created from a quasi random num- 
ber sequence in five dimensions (three coordinates for position 
and two for direction). In our tests a separate random number 
sequence was used for each emission source and for the scat- 
tering process. However, it is still possible that scatterings de- 
stroy the equidistant sampling and could explain the slow con- 
vergence. This was checked by repeating the calculations for 
model Si assuming pure forward scattering. Absorptions are 
calculated explicitly along the photon path so that results are 
not affected by the positions of the scatterings. Fig. [9] shows 
convergence for un-weighted and weighted sampling (meth- 
ods A and B applied). This shows that in normal Monte Carlo 
calculations the convergence rate does finally become ~ 1 /N 
and after some 10 5 photon packages the accuracy reaches that 
of the weighted calculations. When weighted sampling is used 
quasi random numbers give a slightly lower error level but the 
convergence remained at ~ 1/ \N. It may still be possible to 
improve the uniform sampling of scattered flux. This could be 
accomplished by using separate random number sequences in 
different parts of the cloud but there may be other, more prac- 
tical methods. The possibility of ultimately reaching conver- 
gence proportional tol/N makes further studies necessary. 



4.2. Convergence 

Tests showed that simple extrapolation based on previous tem- 
perature values was able to decrease the overall run times by a 
factor of a few. The method requires that temperature values are 
stored for two previous iterations. This extra memory is, how- 
ever, inconsequential for one-dimensional models and is even 
in three-dimensional models only a small fraction of the total 
memory requirements. 

For optically thick clouds and when self-coupling between 
dust temperatures in different parts of the cloud is impor- 
tant AMC-type methods provided equally good or better re- 
sults. The use of the tridiagonal operator is restricted mostly 
to one-dimensional models, but diagonal operator can be ap- 
plied equally well in three-dimensional models. The use of the 
diagonal operator requires storage of photon escape probabil- 
ities, i.e. one number per cell and frequency. This doubles the 
amount of data produced during the simulation of the radiation 
field. If one takes into account all data that are needed in the cal- 
culations (optical depths for absorption and scattering, source 
functions etc.) the percentual increase is not very significant. 
For tridiagonal operator two additional numbers per cell and 
frequency need to be stored since calculations require infor- 
mation on the coupling between neighbouring cells. This could 
possibly be managed even in three-dimensional case. However, 
there the method is useful only if each cell interacts mainly 
with two other cells. This might be the case, e.g., for cylindri- 
cally symmetric models with thin cylinders or in the case of 
a three-dimensional model built of pieces of spherical shells. 
Depending on the geometry (and the available memory) it may 
be useful to include interactions between several neighbouring 
cells. The implementation is not essentially more complicated 
than for the tridiagonal operator. 

The extrapolation method and AMC-methods introduce 
some overhead the importance of which depends on the num- 
ber of simulated photons packages, cells and frequencies. For 
the extrapolation method the simulation part is identical with 
the default method. The time required for the extrapolation is 
small and since extrapolation is done, for example, only ev- 
ery third iteration the increase in run times per iteration is un- 
noticeable. For AMC-calculations with diagonal operator the 
computations differ again very little and no differences in run 
times per iteration was observed. For the tridiagonal operator 
the situation is somewhat different. In our implementation dust 
temperatures were solved in each iteration from a set of non- 
linear equations and this solution was also an iterative one. In 
spite of this it was difficult to see any difference in run times 
per iteration. For example, with model N2 and even in the first 
iterations when the solution was far from the correct one, and 
more sub-iterations were presumably needed to solve the set of 
equations, the difference in the total run time per iteration was 
at most ~1%. 

4.3. Comparison with earlier studies 

In the case of the model Si from Stamatellos & Whitworth 
l2003l and the model Ni from Niccolini et al l2003l our results 
are in perfect agreement with the results given in those articles. 
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Niccolini et al. presented timings for their calculations of the 
model Ni. They found an rms-error below 1% when number 
of photon packages was a few times 10 5 while in our Fig.0]the 
rms-error for the un-weighted Monte Carlo method was still at 
least above 10%. The difference is probably caused by a dif- 
ference in the radial discretization. Niccolini et al. 2003 used 
65 shells and although we have only 51 shells the innermost 
shells have relatively much smaller volumes. Therefore, if nor- 
mal Monte Carlo sampling is used the temperature errors re- 
main very large in these shells and the rms-errors of Figure- 
fleet this fact. When weighted Monte Carlo was used the effect 
of the different volume sizes is removed and 1 % accuracy was 
reached with ~ 10 4 photon packages per iteration. Comparison 
of run times is even more uncertain because of the different 
platforms used. For their method 2 and for 2.4 x 10 6 photons 
Niccolini et al. reported a run time of 500 s on a Cray T3E- 
1200 parallel computer with 16 processors. Since an accuracy 
of 1% was reached with a few times 10 5 photons those compu- 
tations would have taken ~ 50s. With accelerated Monte Carlo 
and weighted sampling our calculations took a total of 12 sec- 
onds on a 600MHz PHI computer. It is not clear whether the 
run times given by Niccolini et al. were for one iteration or the 
whole run. Nevertheless, the comparison seems to indicate that 
the use of weighted sampling and AMC-type acceleration can 
indeed result in significant savings. 

The method of Bjorkman & Wood (2001 ; in the following 
BW) was described to solve the dust temperatures without it- 
erations and therefore the equilibrium temperature calculations 
would require no more time than a pure scattering model. The 
BW method does, indeed, have several advantages over the de- 
fault Monte Carlo method. Firstly, all simulated photon pack- 
ages are used in the derivation of the final solution. In normal 
Monte Carlo runs the radiation field is estimated on each iter- 
ation independently using a new set of photon packages and, 
compared with the total number of simulated packages, the 
random errors are larger. Secondly, in the BW scheme correc- 
tions are applied immediately when a new photon is absorbed. 
In normal Monte Carlo, cells continue to send photons corre- 
sponding to the old temperature until all temperatures are up- 
dated at the end of the iteration. Finally, in the BW scheme the 
photons follow the actual flow of energy so that more emissions 
(and re-emissions) take place in the warmer region. Therefore, 
the sampling of these regions is particularly good. The immedi- 
ate re-emission mechanism helps to ensure that there will be no 
poorly sampled regions that could separate energy sources from 
the rest of the model (see Sect. 13. 1> . When a short- wavelength 
photon is absorbed (e.g., close to a radiation source) this is 
usually followed by a re-emission at a longer wavelength for 
which the optical depths are likely to be considerably lower. 
The resulting increase in the photon free path makes the simu- 
lation procedure much more efficient in the case of large opti- 
cal depths. However, most of the other listed deficiencies of 
the normal Monte Carlo method can be alleviated by using 
weighted sampling and a reference field. 

In the BW method energy conservation is enforced for each 
photon package separately. Once the simulation of the selected 
number of photons from the various radiation sources has been 
completed the dust temperatures will also have reached - apart 



from sampling errors - their final, correct values. It was ar- 
gued that in the case of high optical depths the method does 
not suffer from similar slowdown as A-iteration methods. This 
does not mean that run times would not be affected by r. 
Consider, for example, a cell which is optically thick even for 
the re-emitted radiation. Successive re-emissions and absorp- 
tions slowly increase the temperature of this one cell while the 
flow of information across cell boundaries is minimal. Each 
absorption implies a new calculation of the dust temperature so 
that, like in A-iteration methods, new dust temperatures must 
be solved numerous times. The efficiency of these temperature 
updates is clearly crucial for the overall run times. The BW im- 
plementation used large arrays of pre-tabulated Planck mean 
opacities to speed up the solution (see Bjorkman & Wood l200Tl 
Eq.6). 

In the default Monte Carlo scheme information is first gath- 
ered from a large number of photon packages and new dust 
temperatures are calculated only at the end of an iteration. 
Therefore, the ratio between the number of simulated photon 
packages and the number of temperature updates is much larger 
than in the BW method. This will become significant if calcula- 
tions include several dust populations, grains size distributions 
or transiently heated grains. Against this background, it is in- 
teresting to see a comparison of the methods in the simplest 
case of one dust population, one grain size and grains at one 
equilibrium temperature. We compared our run times with the 
program MC3D 1 by S. Wolf (2003 1 which includes an imple- 
mentation of the BW method. The cloud model consists of cen- 
tral black body source with radius Rg and temperature 2500 K. 
Beyond a central cavity the surrounding cloud extended from 
3i?s to 300 i?s and had a density distribution n ~ r~ 2 . The 
dust consists of 0.12 /im astronomical silicate grains (Draine 
& Lee 1984). We consider two models where optical depth to 
cloud centre is either 100 or 1000. 

In runs with the MC3D program only the number of sim- 
ulated photons was changed and otherwise default parame- 
ters were used (including code acceleration that neglected very 
small fluxes). In our own program we used weighting scheme 
B, a reference field (see Appendix IA.4t and iterations were 
accelerated with the AMC method using a diagonal operator. 
Both the number of iterations and the number of photon pack- 
ages per iteration were changed between runs. The rms errors 
were estimated by comparing obtained dust temperatures with 
reference solutions that were computed separately for both pro- 
grams using a much larger number of photon packages (and it- 
erations). The rms-errors as the function of run times are shown 
in Fig. All runs were made on the same computer (AMD 
Athlon MP 2000+) on a single processor. 

Comparison of the methods themselves is affected by the 
usual uncertainties (programs use different programming lan- 
guages and the final run times depend on details of the actual 
implementation and compiler optimizations). However, some 
conclusions can be drawn. First of all, with the aid of a ref- 
erence field and AMC acceleration the normal Monte Carlo 
method is competitive with the BW method. Secondly, the 

1 The program MC3D can be downloaded at 
http://www.mpia-hd.mpg.de/FRINGE/SOFTWARE/mc3d/ 
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Fig. 10. Comparison of run times between our program (solid 
squares) and the Bjorkman & Wood (2001 1 method (triangles) 
as implemented in the MC3D program of Wolf (2003 1. The 
figures show rms-errors of the computed dust temperatures as 
functions of CPU-time. The two frames correspond to mod- 
els with central optical depth Ty = 100 and Ty = 1000. The 
dotted lines indicate run times of MC3D program when the 
longer initialisation time (pre-calculation of Planck mean opac- 
ities etc.) has been subtracted. 



AMC method seems to have an advantage at higher optical 
depths, in this case at r ~ 1000. As noted above, the BW 
method requires frequent temperature updates. If a normal 
Monte Carlo run is accelerated with an AMC scheme the to- 
tal number of temperature updates remains much lower. By in- 
ference, in situations where temperature updates are more ex- 
pensive the combination of an AMC scheme and the use of 
a reference field may prove to be more efficient than the BW 
method. 

In addition to run times one may compare memory require- 
ments of possible implementations. If, as in the case of the pre- 
vious example, dust properties are identical in all cells the in- 
formation of the general K sca (f) and K a bs(f) curves can be 
kept in memory. In the BW scheme the local density and tem- 



perature would need to be stored for each cell but the r a b s 
and r sca values can be calculated based on the local density 
and the interpolated or calculated n sca (v) and K a b s (^) values. 
Total memory requirement corresponds to only ~ 2 x N ce i\ s 
floating point numbers. Let us next considers a more com- 
plex problem with A pop dust populations, each discretized 
into N s i ze size intervals. In this case memory is required for 
~ A ce ns x iVp p x A s j zes temperature values. If dust abundances 
are not constant, there will be further N cc n s x A pop abundance 
values. Optical depths could be calculated in real time sum- 
ming over the dust populations although this could have some 
negative impact on the run times. However, it seems probable 
that during each temperature update one needs to have access 
to all absorption cross sections so that the absorbed energy (cal- 
culated based on total optical depths) can be divided between 
grains in different grain populations and size intervals. In or- 
der to relax the memory requirements one could read and write 
temperature values directly to disk files but, since these values 
are updated after each absorption, this would cause a signifi- 
cant increase in the run times. Finally, if the problem includes 
transiently heated particles those must be discretized into Ne 
enthalpy bins. No implementation of the BW scheme yet ex- 
ists for non-equilibrium dust grains. However, in that case not 
only are the temperature updates much more time consuming 
(relative to the simulation part) but involve significantly larger 
amount of data, up to JV ce rj s x A pop x A S j ZC x Ne temperature 
values. 

In normal Monte Carlo the simulation of different discrete 
frequencies can be made sequentially and completely indepen- 
dently. Each time the simulation of a new frequency starts a ta- 
ble of the r a b s and r sca values can be calculated for this one fre- 
quency and all cells. A priori, no assumption need to be made 
that absorption and scattering cross sections would be the same 
in all cells. Counters are also needed for photons emitted and 
absorbed at the current frequency. Therefore, the total memory 
requirement is ~ 4 x JV ce rj s numbers. The AMC method with 
diagonal operator adds to this A cc n s numbers and, alternatively, 
a reference field ~ 2x iV ce u B numbers. In runs with both AMC 
(diagonal operator) and a reference field is ~ 7 x Amelia floating 
point numbers i.e. a few times the amount needed in the BW 
scheme. The program version used in this paper follows this 
scheme apart from the fact that absorption an emission coun- 
ters are all kept permanently in main memory (this applied also 
to the additional counters used with the reference field and the 
AMC methods). This avoids the need to read and write tables 
of n a bs and n cm jt values between simulations of different fre- 
quencies. This does not affect the run times of ID models since 
all values fit in disk cache. We tested the effect on run times in 
a more realistic setting, using a 3D model where each cell was 
hit by ~ 100 photon packages per iteration and simulated fre- 
quency. The version where all n a b s and n em it values were per- 
manently in main memory was faster by ^7 %. This shows that 
the overhead associated with these disk operations is not very 
significant. This is not surprising since each value is read and/or 
written once per iteration but is consequently used in the simu- 
lation of possibly hundreds of photon packages. Finally, let us 
consider a case where any of the factors N S i ze , N pop , or Ne 
are larger than one. During the simulation one needs only total 
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absorption and scattering cross sections for one frequency at a 
time and the counters, n a b s , n e mit etc., register the total num- 
ber of events at that frequency. Therefore, memory requirement 
of the simulation are not affected by any of the factors N size , 
7V pop , or N E . 

The traditional Monte Carlo scheme allows an efficient use 
of external files and each value (n cm ;t, n a b s , etc.) is read from 
and written to an external file exactly once per iteration. The to- 
tal overhead from all file operations depends on the number of 
photon packages simulated but should be clearly below 50%. 
In the test mentioned above the overhead was 7% per file of 
iV ce n elements. The reference field and AMC (diagonal opera- 
tor) require, in addition to n a bs and n em it counters, three such 
files. It might also be possible to calculate cumulative absorbed 
energy without storing absorbed energy at each frequency sep- 
arately. This would further eliminate the need for storing n a b s 
values separately for each frequency. In the BW scheme the 
random selection of the frequency of the re-emitted photons 
means that optical depths must be calculated in real time or 
large arrays are needed for pre-calculated values. Each absorp- 
tion event leads to a temperature update and this leads either 
to very frequent disk operations (if temperatures are kept on 
disk) or large memory requirements. In Fig.^|the comparison 
of run times was based on one type of grains at one equilib- 
rium temperature. In more complex problems the ratio of run 
times will probably change (see above) and the ratio between 
the memory requirements of a BW program and our program 
(AMC and reference field included) will probably be of the or- 
der of- (xA pop x N sizc x N E )/3. 

4.4. Monte Carlo calculation at high optical depths 

Large optical depths cause two problems for radiative trans- 
fer calculations. The sampling of the radiation field becomes 
more difficult as the photon free path becomes very short and, 
secondly, iterations converge much more slowly. The first prob- 
lem can be addressed with weighted sampling and the second 
is at least partially remedied with AMC methods. So far we 
have considered spherical models with central optical depths 
up t — 1000. The dust model is now changed so that optical 
depth is r a b s = 10000 at all frequencies and re-emission at 
longer wavelengths no longer provides an easier escape path 
for the radiated energy. 

In an extremely optically thick cell most of the emitted en- 
ergy is absorbed within the same cell and photons packages 
only rarely move to another cell. This means that one must 
create very large number of photon packages before the en- 
ergy flow between cells can be estimated with any accuracy. 
Weighted sampling gives a solution for this particular problem. 
If photon packages are created close to a cell boundary it is 
much more likely that it will eventually cross it i.e. give infor- 
mation about the actual energy flow. For example, one could 
use exponential probability distribution q — e~ T< * s \ where r(s) 
is the optical depth to the closest cell boundary. However, in 
this case the weights of individual photon packages are random 
variables and the total emission from a cell would fluctuate. 
For optimal results these fluctuations should be corrected, for 



example, by adding a few photon packages for which the sum 
of photons corresponds to the difference between the expected 
value and the number of photons so far simulated. 

In the following we use a simpler approach. We assume 
that all photons further than r a b s = 15 from the closest cell 
border are absorbed in the cell. All simulated photon pack- 
ages are started in the remaining volume where distance to the 
closest border is below r a b s = 15. Weighting takes into ac- 
count the fact that packages are started only in a part of the 
cell volume. Same number of photon packages in sent from 
each cell (weighting scheme E) and the number of true pho- 
tons simulated from each cell does not fluctuate. The fraction 
e~ 15 ps 3 • 10~ 7 is sufficiently small so that the true flux from 
deeper in the cell is insignificant compared with the emission 
from regions closer to cell surface, r(s) = 0—15. If optical 
depth of a cell is increased beyond r a b s = 30 the accuracy of 
the sampling will no longer be affected. Without AMC the con- 
vergence would be extremely slow since the ratio between pho- 
tons coming from neighbouring cells is much smaller than the 
number of photons emitted and absorbed within a cell. AMC 
scheme is used to eliminate the local emission-absorption cycle 
and temperature updates depend only on photons that actually 
cross the cell borders. 

Figure ^2 shows some results from our run (solid line). 
From the dust medium 5000 photon packages were simulated 
each iteration. The emission from the central star was divided 
between 1000 iterations using 500 photon packages per itera- 
tion. This was done to improve the use of the reference field 
(see Appendix IA.4> . The total number of iterations was 1500 
and the run time —15 minutes. For reference we plot a di- 
rect numerical solution for a case where energy transfer be- 
tween cells would depend only on cell area and temperature, 
ocijX T^ ust i (dashed line). This would correspond to infinite 
optical depth. The comparison shows that a reasonable solution 
was indeed found in a relatively short time. The small fluctu- 
ations in the computed dust temperatures further indicate that 
the higher optical depths did not cause numerical problems. In 
this case a similar solution could be reached even with regu- 
lar weighting (scheme B) using a minimum of — 20000 pho- 
ton packages per iteration (solid squares). When sampling is 
restricted to regions close to cell boundaries it is possible to in- 
crease the optical depths further. A solution for r = 10 5 model 
was found in same time and with the same number of photon 
packages as for the r = 10 4 model (5000 photon packages 
per iteration and frequency for the dust emission). The result- 
ing dust temperature distribution was (—within noise) identical 
with the ones shown in Fig. II II 

The previous example shows that the Monte Carlo method 
can, in principle, be used even the case of cells with very large 
optical depths. Photon packages are used only to sample the 
energy flux across cell boundaries and there are not necessarily 
any packages that move across even a single cell. Information 
moves across the cloud at the rate of one cell per iteration 
and run times depend directly on the discretization of the op- 
tically thick region. The weighted sampling must be tailored 
according to the cell geometry. In regular geometries (e.g., a 
3D cartesian grid) the implementation is relatively easy and the 
weighting can always be restricted to regions where it is actu- 
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Fig. 11. Dust temperatures in a model with a central black 
body source and a surrounding dust cloud with an optical depth 
r = 10000 at all frequencies. The solid line shows the results 
when emission was simulated only from those regions where 
distance to the closest cell boundary was below r a b s = 15 
(see text). The solid squares correspond to a run with normal 
weighting (scheme E) and a factor of 4 larger number of pho- 
ton packages. Direct numerical solution for a case where en- 
ergy transfer across a surface is oc A x T$ ust (infinite optical 
depth) is drawn with a dashed line. 

ally needed. In optimal case the optically thick region would be 
treated separately so that, in addition to weighting, even com- 
plete iterations could be carried out independently from the rest 
of the model. Such methods are already in preparation for hier- 
archical models. Hierarchical discretization makes it possible 
to keep optical depths of individual cell lower and, if incom- 
ing intensity is saved at grid boundaries, the calculations of the 
various sub-grids can be carried out relatively independently. 
Consequently, a large number of iterations required in certain 
sub-grids does no longer necessarily have a large impact on the 
overall run times. 

5. Conclusions 

We have discussed various methods that can be used to improve 
the efficiency of Monte Carlo radiative transfer calculations for 
dust scattering and emission. The accuracy of the simulation 
of the radiation field can be improved with weighted sampling 
and/or by employing a reference field. To our knowledge this 
is the first time that the reference field and the weighted angu- 
lar distribution of scattered and dust emitted photons are con- 
sidered in connection with dust continuum calculations. More 
surprisingly, even spatial weighting of initial positions of pho- 
ton packages is usually not used. Tests with internally and ex- 
ternally illuminated/heated spherical model clouds lead to the 
following conclusions: 

- improper sampling can lead to wrong results although cal- 
culations have some appearance of convergence 



- weighted sampling resulted in significant savings in run 
times ranging from a factor of ~ 1 — 10 in externally heated 
models to several orders of magnitude in internally heated 
and more optically thick models 

- one should use weighting schemes that reduce fluctuations 
in the sampling, e.g., simulating the same number of photon 
packages from each cell - this reduces sampling errors and 
can lead to a faster convergence, ~ 1 /n p , as function of the 
number of photon packages 

- weighting can be applied even to the angular distribution 
of scattered photons: it can be used to improve sampling in 
selected regions inside the cloud and to lower the noise in 
the scattered flux observed in the selected general direction 

- weighting schemes are easy to implement and do not signif- 
icantly affect the run times per photon package: they should 
be used in most calculations 

- use of quasi random numbers (instead of pseudo random 
numbers) brought some improvement, but only in the case 
of pure forward scattering did the the convergence ap- 
proach l/n p 

- the use of a reference field improves efficiency of calcu- 
lations when several iterations are needed: emission from 
constant radiation sources (e.g., background) needs to be 
simulated only once and the number of other photon pack- 
ages can be reduced leading to very significant savings in 
run times 

The convergence (as function of the number of iterations) 
can be improved by simple extrapolation that is based on dust 
temperature from previous iterations or by using accelerated 
Monte Carlo (AMC) methods. These methods were, for the 
first time, used in Monte Carlo calculations of dust emission. 
Our tests showed that: 

- heuristic extrapolation of temperature values requires some 
extra storage but involves very little computations and can 
reduce run times by a factor of a few 

- accelerated Monte Carlo method where internal emission- 
absorption- cycle is eliminated is equally effective and also 
incurs very little computational overhead 

- for one-dimensional models one should consider including 
interactions between neighbouring cells; the run time per 
iteration remains essentially unchanged but the number of 
iterations is significantly lower for optically thick models 

- heuristic extrapolation can be successfully combined with 
accelerated Monte Carlo methods 

Considering the storage and computational overhead all 
methods (weighted sampling, reference field, heuristic extrap- 
olation and accelerated Monte Carlo) are also suitable for the 
use with two- and three-dimensional models. Tests with ID 
cloud models showed that when weighted sampling and AMC 
methods are used the run times are very similar to those of 
the Bjorkman & Wood J2001t simulation scheme. However, 
for larger and more complex models (e.g., with several dust 
species) the memory requirements are lower and the the run 
times may compare favourably with the Bjorkman & Wood 
method. 
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Appendix A: Implementation of the program 

A.1. The basic program 

The methods discussed in this paper are implemented in the ra- 
diative transfer program used in Juvela & Padoan (2003 1. The 
simulated spectral range is covered by a discrete set of frequen- 
cies. At each frequency the radiation field and the resulting dust 
absorptions are simulated with a number of photon packages. 
These are sent separately from discrete sources (e.g., the central 
star in models Ni and N2), from the isotropic background and 
from the dust filled cloud volume. In the basic method the true 
number of photons emitted from each of these sources is di- 
vided equally between photon packages sent from that source. 

On the stellar surface and on the border of the spherically 
symmetric cloud the photon packages are created at locations 
{rcos4>isin0i,rsin(j)isin9i,rcos6i), where r is the radius of 
the object. The random angles 9t and tfii are obtained from 
equations 



acos{u) 1 4>i = 2nu, 



(A.l) 



where u denotes a uniform random numbers in the interval 
[0, 1 [. For dust emission within the cloud volume the radius r is 
replaced with a random variable r.i — M 1 ^ 3 i? c ioud ( see below). 

For emission from stellar surface and the background the 
photon packages have initial directions over 2ir solid angle that 
are determined by random angles 



acos(\/u), <f>i = 2iru. 



(A.2) 



Here 6i is an angle from the normal of the surface and <fti ro- 
tation around the normal. For dust emission the packages are 
uniformly distributed over An radians and angles 9 (measured 
from arbitrary reference direction) are 



acos(u) 



(A3) 



Consider a step for which the optical depth of scattering 
is T scat . After this step the fraction of unscattered photons is 
1 — e~ Tscat . This is also the cumulative probability distribution 
for the photon free paths and the inverse functions gives the 
formula for simulation of random free paths, 



-ln(u) 



(A.4) 



In this equation u is again a uniformly distributed random num- 
ber. When photon package moves through a cell the optical 
depth of absorption, r a b s , is calculated and the number of ab- 
sorbed photons is removed from the photon package and added 
in an absorption counter in the cell. If photon package has ini- 
tially photons the number of absorptions is 

An v = <(1 - e- T °- b *). (A.5) 

When the total optical depth of scattering reaches r° CQ the pho- 
ton package is scattered toward a new direction that is deter- 
mined by the dust scattering function. A new value of t° cq is 



generated and simulation continues until the number of photons 
remaining in the photon packages has become insignificant or 
the package exits the cloud. If package did scatter before exit- 
ing it can be registered as a sample of the scattered flux. 

In an alternative simulation scheme the free path is calcu- 
lated as above, t° = —ln(u), and scattering occurs when total 
optical depth along the path reaches this value. Absorptions are 
not calculated along the photon path. When package scatters 
the fraction of photons corresponding to the dust albedo is re- 
tained in the package and the rest are added to the absorption 
counters of the current cell. This scheme is faster when opti- 
cal depth is large but leads to poor sampling in regions of low 
optical depth. This method was not used in this paper. 

In some runs the method of forced first scattering (e.g., 
Mattila fT970l was also used. In that method one first calculates 
the optical depth ET sca to the cloud edge along the original di- 
rection. The fraction e~ Tsca of moves along this line without 
scattering. The remaining fraction, 1 — e~ Tsca , does scatter at 
least once and a conditional probability of their free path will 
be calculated: 'where will a photon scatter if it does scatter be- 
fore the cloud border' . The normalized cumulative probability 
density function is 



P(r) = 



1 - e 



(A.6) 



1 _ e -Sr 8ca ' 

and random free paths are calculated using function P _1 , 

**•«)). (A.7) 



r» ca = -ln{l-u{l-e 



The scattered part is calculated explicitly for each package and 
this improves accuracy of the estimated scattered flux, espe- 
cially if the optical depth is low. 

A.2. Weighted sampling 

In the following function p(r, ft) denotes the original proba- 
bility density function that in normal Monte Carlo runs de- 
termines the number of photon packages sent from different 
position, r, and towards different directions, ft. The photon 
packages can be re-distributed according to another probabil- 
ity function q(r, SI). The probability distribution of the actual 
photons must not be changed and this is accomplished by mul- 
tiplying the number of true photons included in the package by 
the ratio 



W 



p(r,Q) 



(A.8) 



If number of packages is increased in some region (q > p) the 
weight of that package is correspondingly decreased (W < 1). 



A.2.1. Method A 
packages 



positions of external photon 



For external photons entering a spherical cloud the original 
probability distribution of the impact parameter d is p(d) = 
2(d/i? c ioud), where i? c ioud is the cloud radius. The cumula- 
tive probability distribution is P(d) — (d/R c i ou d) 2 and ran- 



dom impact parameters are di 



P _1 (u) = sfuR, 



loudi 
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where u is a uniform random number. In addition to d one 
needs a random angle 4> € [0, 2ir] to determine a rotation 
in the plane perpendicular to the direction of the incoming 
photon package. Uniform sampling is not adequate if the in- 
nermost shells are very small. The sampling of the inner re- 
gions could be improved, for example, by using a centrally 
peaked probability distribution q(d) = (1 + a)(d/i? c ioud) a 
with a < 1. In this case the cumulative probability density 
function is Q(d) — {d/ R c \ ou d) 1+a and the inverse function 



gives formula d, = R c \ ou dU 



(l+a)- 



(a ^ —1) for generating 



impact parameters. Usually one would include in the package 
n v actual photons corresponding to the total number of pho- 
tons entering the cloud divided by the number of photon pack- 
ages generated from the background. Now this number must be 
scaled with W = p(d i )/q(d l ), 



P{dj) 

q(di) 



i 



pa— 1 jl- 
~- rl cloud"i 



(A.9) 



In this paper we use a scheme where identical number of pho- 
ton packages is sent towards each annulus as defined by the ra- 
dial discretization. The annulus is selected systematically and 
within the annulus the probability follows the default p(d) ~ d 
dependence, 



di 



Rk- 



R 



loud 



•«[( 



Rk 



loud 



Rk-i 
Rc 



loud 



(A. 10) 



The index k refers to the selected annulus and Rk is outer ra- 
dius of the shell k. The photon numbers are scaled with the 
ratio between the probability of selecting this annulus in the 
normal case and in the new scheme. The first probability is 
proportional to the area of the annulus and, since the annulus is 
selected systematically, the latter is simply one over the number 
of shells, 1/iVghoiis- The weight is 



Wi = [( 



Rk 



R, 



cloud 



■) 2 -( 



R k - 



R 



cloud 



-) 2 ] x N t 



hells- 



(A.ll) 



A.2.2. Method B 
packages 



positions of internal photon 



For emission within the cloud volume the original cumulative 
probability distribution for distances from the cloud centre is 
P{r) — (?"/i?cioud) 3 and usually distances would be gener- 
ated from formula i\ = P _1 (?i) = i? c ioud it 1 / 3 . We are using 
a scheme where shells are selected systematically and within 
each shell we use the normal probability distribution, ~ r 3 . 
After selection of the shell k the actual radius for the emission 
event is generated from equation 



(3hzL)3 + u [(_?*_)3 _ (^*=L)3) 



R 



cloud 



R 



cloud 



R, 



(A. 12) 



cloud 



The weight of the photon package is equal to the the cell vol- 
ume (relative to cloud volume) multiplied by the number of 
shells, 



[( J^_)3 _ (JhzLf] x N, 



R, 



cloud 



R, 



shells • 



(A.13) 



cloud 



A. 2. 3. Methods C and D - angular distribution 

Weighting can be applied to angular distributions of emitted 
and scattered photons. For emitted photons the original distri- 
bution is uniform while for the scattered photons the original 
distribution p(9, fa) is determined by the scattering function. 
The scattering function is calculated in a coordinate system 
fixed by the original direction of the package. In both cases 
we re-distribute photon packages according to an distribution 
q(9) oc e~ lb . In this equations 9 is an angle from a selected 
direction T which can point, e.g., towards the observer or the 
cloud centre. After normalization this function becomes 



9(0) 



1 - e-T* 

and the cumulative probability density function is 



Q(0) 



l 



-70 



1 - e-~i* 



(A. 14) 



(A.15) 



Function Q 1 is used to generate random angles from distribu- 
tion q{9), 



= — ln[l-u(l- 



(A.16) 



The rotation around the direction toward T is specified by a 
uniformly distributed random angle <\>i € [0, 2w]. In normal 
Monte Carlo simulation the distribution of emitted packages 
is isotropic, p(9) — ^sin9, and the weight is obtained from 
Eq. IA.8l as p(9)/q(9). For scattering the original angular distri- 
bution is not uniform and we have weights 



Wi = 



(A.17) 



The scattering function (i.e. function p) involves angles and 
4>i that are defined relative to the original direction of the pack- 
age while in the function q the angles 9i and fa are relative to 
the selected direction T. Functions p and q are uniform with 
respect to the angles <f>' and fa respectively. However, because 
of the rotation between the two coordinate systems the value of 
9' depends on the selected angle fa and the weight W depends 
the selected values of both 9 and (j>. 

A.3. Accelerated Monte Carlo 

The basic idea behind Accelerated Monte Carlo (AMC) meth- 
ods is that part of the radiative interactions are treated explicitly 
when dust temperatures are updated. The methods are com- 
pletely analogous with ALI methods but implementation dif- 
fers since we compute numbers of absorbed photons and not 
intensities. 

Each cell has counters for the number of photons that are 
absorbed in that cell during one simulation of the radiation 
field. New dust temperature estimates are calculated by balanc- 
ing the absorbed energy with emission that depends on the dust 
temperature. This calculation takes into account the change in 
the emitted energy but not the effect that a change in the tem- 
perature has on absorptions. In an optically thick cell the ab- 
sorption counters contain mostly photons that were originally 
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emitted within the same cell. If temperature is, for example, in- 
creased the calculations assume that all added emission escape 
the cell. In reality most photons are re-absorbed and the net 
flow of energy out from the cell is much smaller. Consequently, 
the needed temperature correction is severely underestimated 
and convergence remains slow. 

In AMC this problem is removed by treating explicitly 
those photons that are absorbed within the same cell from 
which they were emitted. The simulation proceeds in normal 
fashion except that absorbed photons are counted separately 
depending on whether they were originally emitted from the 



,abs 



same cell or not. We denote these counters with rtf nt s and n^ xt 
the indices referring to the origin of the photons (internal vs. 
external). In normal runs the absorption counters (number of 
photons per unit frequency interval) correspond to the sum of 
these two counters and the dust temperatures is solved from 
an equilibrium condition for unit volume, 



n abs (v)hvdv = 4tt / K v B v {T A )dv 



(A.18) 



In AMC schemes this equation will be modified. In addition to 
n fnt S we nave a ^ so information about the number photons emit- 
ted from the cell, n emlt , and we can calculate photon escape 
probabilities 



= 



^ emit abs 



(A.19) 



On the right hand side the part 1 — represents photon emitted 
and absorbed in the same cell. On the left side this corresponds 



to counters nf^. This part can be subtracted and we are left 
with a condition 



abs 



hvdv — An I K, v B v (T&)0dv. 



(A.21) 



The left hand side is the net inflow of energy into the cell, the 
right hand side is the net energy loss and the cycle of local ab- 
sorptions and emissions has been eliminated. This form will 
lead to faster convergence, especially if escape probability (3 
was small. In ALI terms this corresponds to a separation of 
the diagonal part of the lambda operator. The lambda opera- 
tor is defined by the relation J = AS and diagonal elements 
describe the effect of the local source function on the local 
mean intensity. In AMC this separation was accomplished by 
using separate counters, nf^. In our implementation the tem- 
perature is solved from Eq. IA.21I with Newton-Rhapson iter- 
ation. In large models the run times could be decreased by 
pre-calculating a direct mapping between incoming energy (left 
side of Eq. IA.2H and dust temperature. That way the solution 
of Ea. lA.2ll would involve only a single table look-up (and pos- 
sibly an interpolation). However, Ea. IA.2ll represents a small 
fraction of all computations and the effect on run times would 
be negligible. 

The AMC methods can be extended by considering radia- 
tive interactions between cells. Additional counters can be used 



to register separately those of the absorbed photons that were 
emitted from an immediate neighbour. In a spherical model this 
means the neighbouring shells with indices k — 1 and k + 1. 
The total number of absorbed photons can be written 



abs 



abs 



abs 



abs 



^abs 



n k, int + n k,cxt + n k-\M + n fe+l,fc; 



(A.22) 



where terms nfj s correspond to photons emitted in cell i and 
absorbed in cell j. This time the counter includes only 
photons originating beyond the neighbouring cells. Together 
with the number of emitted photons these define ratios 



, abs 



(A.23) 



for i = j ± 1. One can also define £ k .k which is equal to 1 — 0. 
The terms ^ j describe the radiative interaction between cells: 
what part of photons emitted in cell i is absorbed in cell j. 
Actually, & j is based on photon numbers per unit volumes so 
that in the following the cell volumes are not explicitly vis- 
ible. In the equilibrium equation counters k can be re- 
placed with calculated emission from shell k ± 1 multiplied 
with the corresponding factor £k±i,k- The escape probabilities 
are again used to eliminate local absorptions, and the equi- 
librium condition for cell k becomes 



The equilibrium equation can be re-written in the form 

J (n^ t s +n^t) hvdv = AnJ K v B v (T d )[(l-0)+0]dv.(A.2O) + Jj 



+4tt / K k -x{v)B v {T*- l )^i >k {v)dv 
- : -4;r / ,; k+1 (v)B v (T^ +1 )^ k+1:k (v)dv 
An f K k {v)B v (T%) 0(v)dv 



(A.24) 



The right hand side identical with Eq. IA.2ll The left hand side 
depends now explicitly on the dust temperatures of the neigh- 
bouring cells and all temperatures must be solved from a non- 
linear set of equations. This has the advantage that tempera- 
ture changes in neighbouring cells are taken into account im- 
mediately and not only after next simulation step. If shells are 
optically thick there will be little radiation coming from be- 
yond the neighbouring shells (nf 3 ^ w 0) and the solution of 
the Ea. lA.24l is immediately close to the true solution. In ALI 
terms this corresponds (in one-dimensional case) to the use of a 
tri-diagonal operator. In our implementation Newton iterations 
are used to solve for one temperature while other temperature 
values were kept constant. This is repeated for all cells until 
solution is fully converged. The procedure is not particularly 
efficient but still fast enough so that overall run times per itera- 
tion were practically the same as when no AMC methods were 
used. 

Since coupling between all cells is not stored the simula- 
tion of the radiation field must be repeated for each iteration 
in the usual fashion. The photon escape probabilities and the 
coupling between neighbours remain, in principle, unchanged. 
We repeat, however, also these steps so that the calculations re- 
main consistent if a different set of random numbers is used on 
different iterations. 
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Further AMC methods can be constructed by consider- 
ing interactions between more cells. For example, in a three- 
dimensional cartesian grid one might consider interactions be- 
tween a cell and its six neighbours along coordinate axes. 
However, each interaction requires separate counters n abs and 
the total memory requirements of the program would be in- 
creased by a factor of a few. 



A.4. Reference field 

Normally Monte Carlo sampling is used separately on each it- 
eration to estimate the strength of the total radiation field. If 
emissions and absorptions are known for some reference situ- 
ation, we need to estimate only the difference between the true 
and the reference field and sampling errors are correspondingly 
smaller. Normal Monte Carlo calculations are based on the cur- 
rent estimate of the number of photons emitted from each cell, 
n emlt , and the number of photons absorbed in each cell, n abs , 
which results from the simulation of the radiation field. We use 
as the reference the solution from the previous iteration. Before 
dust temperatures are updated the numbers are are copied as 
counters for the reference field 



-,abs 
*ref 



^abs 



'ref 



As new dust temperatures are solved we get new estimates for 
n emit an( j tne f u ow i n g simulation is used to determine the 
change from the current reference situation i.e. from the pre- 
vious iteration. Therefore, simulated photon packages corre- 
spond to differences 

An omit = n omit - n^fS (A.26) 

and the resulting absorptions, An abs , constitute a correction 
to the number of absorptions found in the reference case. The 
closer the reference field is to the actual field the smaller are 



the At 



and A , 



,abs 



terms and the smaller the noise in the 



estimated difference between the 'true' field (corresponding to 
the current temperature estimates) and the reference field. Note 
that An cmlt and An abs can also be negative. After the sim- 
ulation new temperatures are calculated using the sum of the 
reference field and the corrections 



^abs 



abs 
Vef 



An abs 



(A.27) 



and at the same time, before new values for n cmlt are calcu- 
lated, the description of the reference field is updated 



n abs 
*ref 



abs 
l rcl 



An 



abs 



emit 
A ref 



_ emit 



(A.28) 



As the reference get closer to the true solution also the fluctua- 
tions from the simulation of the difference field become smaller 
and, in principle, the final noise depends on the total number of 
photon packages and not just on the number of packages sim- 
ulated during one iteration. Therefore, given noise level can be 
obtained using a lower number of photon packages per itera- 
tion. Another advantage is that emission from background and 
discrete sources need to be simulated only once. On the follow- 
ing iterations their effect is already included in the reference 
field. 



The previous scheme is not optimal in the sense that iter- 
ations before large temperature corrections have larger impact 
on the final solution. In an extreme case, if solution converges 
after one iteration the current fluctuations are 'frozen' in the 
solution which does not change on any of the following itera- 
tions. This problem can be avoided by using as the reference 
a running average of the n abs and n cmlt where early iterations 
(when the solution is still far from the correct one) are given 
a smaller weight. Alternatively, one can divide the emission 
of the photons from heating sources over many iterations thus 
avoiding a large temperature jump on early iterations. As the 
solution converges at a nearly uniform pace all iterations con- 
tribute equally to the estimated reference field and the final 
noise level is correspondingly lower. In practise, this method 
was found to work very well. It does slow down the conver- 
gence but the impact is probably no more than ~50%. The 
method was used in Sect. 14.31 where emission of the photons 
from the central source was divided equally over ^2/3 of the 
total number of iterations. 

A reference field can be used together with AMC methods. 
As long as optical depths are temperature independent each it- 
eration gives additional samples for the estimation of escape 
probabilities and the strength of the radiative coupling between 



(A.25) cells. Equations lA27tiA28l can be applied to counters 



.abs 



,abs 



etc. or one can directly calculate a suitably weighted run- 
ning averages of the (3 and £jj factors. 
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